Optical analogue of relativistic Dirac solitons in binary waveguide arrays 
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We study analytically and numerically an optical analogue of Dirac solitons in binary waveguide 
arrays in presence of Kerr nonlinearity. Pseudo-relativistic soliton solutions of the coupled-mode 
equations describing dynamics in the array are analytically derived. We demonstrate that with the 
found soliton solutions, the coupled mode equations can be converted into the nonlinear relativistic 
ID Dirac equation. This paves the way for using binary waveguide arrays as a classical simulator 
of quantum nonlinear effects arising from the Dirac equation, something that is thought to be 
impossible to achieve in conventional (i.e. linear) quantum field theory. 



Introduction - Waveguide arrays have been used 
intensively to simulate the evolution of a nonrelativis- 
tic quantum mechanical particle in a periodic potential. 
Many fundamental phenomena in nonrelativistic classi- 
cal and quantum mechanics such as Bloch oscillations 
[U [2] , Zener tunneling [3J H] , optical dynamical localiza- 
tion [5], and Anderson localization in disordered lattices 
[H] have been simulated both theoretically and experi- 
mentally with waveguide arrays. In a recent study it was 
shown that, rather surprisingly, most of nonlinear fiber 
optics features (such as resonant radiation and soliton 
self-wavenumber shift) can also take place in specially 
excited arrays [7J. Recently, binary waveguide arrays 
(BWAs) have also been used to mimic relativistic phe- 
nomena typical of quantum field theory, such as Klein 
tunneling 8 , 9J , the Zitterbewegung (trembling motion of 
a free Dirac electron) [101 111) , and fermion pair produc- 
tion [TJ], which are all based on the properties of the 
Dirac equation [T3]. Although there is as yet no evi- 
dence for fundamental quantum nonlinearities, nonlinear 
versions of the Dirac equation have been studied since 
a long time. One of the earlier extensions was provided 
by Heisenberg [T3J hi the context of field theory and was 
motivated by the question of mass. In the quantum me- 
chanical context, nonlinear Dirac equations have been 
used as effective theories in atomic, nuclear and gravi- 
tational physics [TSTfTS] and, more recently, in the study 
of ultracold atoms [TJO HO]- To this regard, BWAs can 
offer a rather unique model system to simulate nonlinear 
extensions of the Dirac equation when probed at high 
light intensities. The discrete gap solitons in BWAs in 
the classical context have been investigated both numer- 
ically [2D - f2"5] and experimentally [21] • In particular, in 
Ref. [22 soliton profiles with even and odd symmetry 
were numerically calculated and a scheme with two Gaus- 
sian beams, which are tuned to the Bragg angle with op- 
posite inclinations, was proposed to efficiently generate 
gap solitons. In Ref. [53] solitons were experimentally 
observed when the inclination angle of an input beam is 



slightly above the Bragg angle. 

Inspired by the importance of BWAs as a classical sim- 
ulator for relativistic quantum phenomena, and also by 
past achievements in the investigation of discrete gap soli- 
tons in BWAs, in this Letter we present analytical soliton 
solutions of the discrete coupled-mode equations (CMEs) 
for a BWA and construct Dirac solitons of a nonlinear rel- 
ativistic ID Dirac equation in the quasicontinuous limit. 
This paves the way for using BWAs to simulate nonlin- 
ear extensions of the Dirac equation that violate Lorentz 
invariance [35] , as well as other solitonic and nonsolitonic 
effects of nonlinear Dirac equations. 

Analytical soliton solutions — Light propagation in a 
discrete, periodic binary array of Kerr nonlinear waveg- 
uides can be described, in the continuous-wave regime, 
by the following dimensionless CMEs [H [21] : 
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where a n is the electric field amplitude in the nth waveg- 
uide, z is the longitudinal spatial coordinate, 2a and k 
are the propagation mismatch and the coupling coeffi- 
cient between two adjacent waveguides of the array, re- 
spectively, and 7 is the nonlinear coefficient of waveg- 
uides which is positive for self-focusing, but negative for 
sclf-dcfocusing media. For simplicity, here we suppose all 
waveguides have the same nonlinear coefficient, but even 
if these nonlinear coefficients are different (provided they 
are comparable), then analytical soliton solutions shown 
later will not be changed, because as explained later, one 
component of solitons is much weaker than both unity 
and other component, and thus one can eliminate the 
nonlinear term associated with this weak soliton com- 
ponent. In the dimensionless form, in general, one can 
normalize variables in the above equation such that 7 and 
k are equal to unity. However, throughout this Letter we 
will keep these parameters explicitly in Eqs. ([I]). Before 
proceeding further, it is helpful to analyze the general 
properties of the general solutions of Eqs. (01). First of 
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all, let us assume that (a2 n ,a 2n _i) T 



i 2n (<P2n,<P2n-l) T 



is one solution of Eqs. with if2n an d <f2 n -i being ap- 
propriate functions. In this case, if we change the sign of 
7 while keeping the other two parameters constant, one 
can easily show that a new solution of Eqs. will be 
(a2n,ci2n-i) T = i 2n (<Pz n - 1 ' ^*2n) T ^ where * denotes the 
complex conjugation. Secondly, if the sign of a is changed 
while other parameters kept constant, then a new solu- 
tion of Eqs. will be (a 2n , a 2 n-i) T = « 2 ™(^2n-l, </?2n) T - 
Of course, when a changes sign, we still have the same 
physical system, but with a shift of the wavenumber po- 
sition n in Eqs. by one. The above simple rules allow 
us to quickly find other solutions and their symmetries if 
one particular solution is known, as will be shown later. 

In the specific case when all three parameters 7, k, and 
a are all kept positive, we look for analytical solutions of 
motionless solitons of Eqs. in the following form: 
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where n £ M. characterizes the beam width (i.e. the 
average number of waveguides on which the beam ex- 
tends), and coefficients 6, d and / are still unknown. In 
the system without any loss or gain of energy (i.e., when 
k, a and 7 are all real), the coefficient / must also be 
real, but b and d can be complex in general. Inserting 
the ansatz into Eqs. , assuming a priori that the 
component d2 n -i is much weaker than both unity and 
the other component a2«, such that one can eliminate 
the nonlinear term for a 2n -\, and also assuming that the 
quasicontinuous limit is valid (i.e. n is large enough), 
after some lengthy algebra one gets: 



fd 
inb 
fb 



ttbi — ad, 
2l\d\ 2 d/nl 
ab + Adni/riQ. 



(3) 
(4) 
(5) 



Extracting / and b from Eq. and Eq. 0, respec- 
tively, then inserting them into Eq. we will get one 
quadratic equation for d 2 , and thus can find the values 
for b, d and /. Note that one needs to keep only solutions 
which satisfy the above assumption that |<22n-i| ^ \o-2n\- 
The final solution in the case when 7, k, a > is: 



a2n(z) 
a2n-l(z) 



.sech(^)e* Z( "o^ 



sech(^i)tanh(^=i)e 

^ Tin / V Tin ' 



iz(- 



-a) 



(6) 

It is worth mentioning that the analytical soliton solu- 
tion in the form of Eqs. is derived under two condi- 
tions: (i) the beam must be large enough such that one 
can operate in the quasicontinuous limit instead of the 
discrete one; and (ii) n \cr\ ^> 2k. The latter condition is 
easily satisfied if (i) is held true and if a is comparable 
to k [TT]. If condition (ii) is not valid, one can still easily 



get the analytical solution for b, d and / from Eqs. 
- ([5]), but they are a bit cumbersome and for brevity we 
do not show it here. The solution in form of Eqs. 
represents a one-parameter family of discrete solitons in 
BWAs where the beam width parameter rin can be arbi- 
trary, provided that no > 4, a surprisingly small number 
for the quasi continuous approximation to be valid. 
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FIG. 1: (Color online) Discrete soliton profiles (a,b) for even 
and odd symmetry, respectively. Full circles mark the field 
amplitudes across the BWA. Parameters in (a): k =1; 7 =1; 
a — 1.2; and no = 5. After getting the even symmetry profile 
in (a), we construct the odd profile in (b) by switching the sign 
of a and following the symmetry transformations explained in 
the text. 

In Fig. [lja) we plot the soliton profile with even sym- 
metry calculated by using Eqs. ([6} at z = with full 
circles marking the field amplitudes across BWAs, for the 
parameters given in the caption. Note that soliton profile 
in Fig. [lja) consists of two components: one strong com- 
ponent ai n and another much weaker component ai n -\ 
[see also Fig. (2jc)]. Once we get the soliton solution in 
Fig. [lja), we can construct another soliton solution of 
the same physical system by changing the sign of a and 
following the rules explained in the previous section. In 
that way we obtain the odd symmetry soliton profile de- 
picted in Fig. [ijb). It is important to mention that in 
the case of self- focusing media (7 > 0) , for both even and 
odd symmetries the strong component is always located 
at waveguides with larger propagation constant [channels 
with +|(t| in Eqs. 0], whereas the weak component is 
located at waveguides with smaller propagation constant 
[channels with — | a\ in Eqs. Q]. We are also able to con- 
struct the soliton solutions for the self-defocusing media 
which also possess soliton solutions with even and odd 
symmetries. The only difference with the self-focusing 
media is that now the strong (weak) component is lo- 
cated at waveguides with smaller (larger) propagation 
constant. 

Soliton propagation and generation - Equation 
and the associated solutions obtained by the above sym- 
metry transformations provide the analytical forms of the 
two discrete gap soliton branches numerically found in 
|22) . We note that the propagation constant / of the 
two solitons, given by / = — a + 2k 2 / \n^a) , falls in the 
minigap of the superlattice, near the edge of the lower 
miniband (because 2k 2 / '(jIqCt) <C cr), and thus they are 
expected to be stable [22] • As an example, in Fig. [2]ja) 
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FIG. 2: (Color online) (a,b) Soliton propagation in the (n, z)- 
plane (a) and its Fourier transform in the (k, z)-plane (b) with 
even symmetry profile at the input, (c) Absolute values of the 
field amplitudes for intense (|d2n| with solid line and square 
markers) and weak (|a2 n -i| with dashed-dotted curves and 
round markers) components of soliton at four different values 
of 2 = (red curves); 50 (blue curves); 140 (green curves); 
and 200 (black curves). Soliton profile is so well preserved 
that all these curves just stay on top of each other and one 
can see only the output black curves, (d) Phase pattern S/tt 
of soliton profiles at four above values of z. Colors of curves 
in (d) have the same meaning as in (c). Parameters: k =1; 
7 =1; a = 1.2; and no = 5. All contour plots are shown in 
logarithmic scale. 



we show the soliton propagation along z as obtained by 
numerically solving Eqs. ([I]) with an input soliton taken 
from Eqs. (j6|) at z = 0, demonstrating that the soliton 
profile is well preserved during propagation. Parameters 
used for Fig. [2] are the same as in Fig. [jja). The evolu- 
tion of the Fourier transform of the field a n in Fig. [2^a) 
along z is shown in Fig. [2jb) where the wavenumber k 
represents the phase difference between adjacent waveg- 
uides. Due to the periodic nature of BWAs, within the 
coupled mode approximation, it suffices to investigate k 
in the first Brillouin zone — tt < k < 7r |26j . One very im- 
portant feature of the wavenumber evolution in Fig. [2jb) 
is the fact that there are two components of wavenumber 
centered at k = ±7r/2 which correspond to two Bragg an- 
gles [TT] with opposite inclinations. These two wavenum- 
ber components are generated at the input and preserve 
their shapes during propagation along z. This feature of 
k indicates that the soliton operates in the region where 
CMEs could potentially be converted into the relativistic 
Dirac equations describing the evolution of a freely mov- 
ing relativistic particle [TUHH] . We will come back to this 
important point again later. Figure [2jc) shows the two 
components of the soliton profile at odd and even waveg- 
uide position n. The strong component with solid curves 
and square markers represents the field profile |a 2n | at 
even waveguide positions, whereas the weak component 
with dashed-dotted curves and round markers represents 



the field profile |d2 n -i| at odd waveguide positions. Field 
profiles in Fig. [2^c) are taken at four values of propaga- 
tion distance z = (red curves); 50 (blue curves); 140 
(green curves); and 200 (black curves) - only the black 
curves are actually visible since the the profile is perfectly 
preserved during propagation with a very high precision. 
The soliton profile also perfectly preserves its phase pat- 
tern across the array [Fig. IJd)]. From Eqs. one 
can easily see that as the waveguide position variable n 
runs, the phase pattern of the soliton must be periodic 
as follows: S n — ...(p, p), (p + 7T, p + vr), (p, p)... where p 
also changes with z. This pattern is only broken at the 
soliton center point where the function tanh in Eqs. ^ 
changes its sign. This phase pattern is shown in Fig. 
[2jd) where different colors with meanings as in Fig. [2jc) 
depict pattern at different values of z. The sequence in 
the phase is important because it allows us to convert 
Eqs. ([I]) into the nonlinear Dirac equation as we shall 
show shortly. Note that the soliton whose propagation is 
shown in Fig. [2] is the one with even symmetry in Fig. 
[T] (a). Our simulations similarly show that the profile 
of soliton with odd symmetry in Fig. [T] (b) is also well 
preserved during propagation, and we have checked that 
this is true even in presence of quite a strong numerical 
noise, demonstrating the robustness and the stability of 
our solutions. 

Although the soliton solutions given by Eqs. (|6j) are 
exact, it is important to consider the possibility to gener- 
ate the new gap solitons by an input beam with a simpler 
(and more experimentally accessible) profile. Due to the 
wavenumber structure shown in Fig. [2]^b), one can inter- 
pret the soliton as a combination of two beams launched 
under two Bragg angles with opposite tilts k — ±tt/2, 
similarly to what was suggested in Ref. [35] ■ Here we 
propose to generate the soliton by an input with a simple 
phase pattern where the phase difference between adja- 
cent waveguides is equal to tt/2 across the array. The 
input condition is taken to be A n = a n exp(imr/2) where 
a n is given by Eqs. jfj} at £ = 0, but without the term 
i 2n . Note that, since |a 2n _i| <C |a2n|> this input condi- 
tion can be approximately achieved by exciting the BWA 
with a broad beam tilted at the Bragg angle, with the odd 
waveguides in the structure being realized at some spatial 
delay Az inside the sample (so as they are not excited at 
the input plane); see the scheme shown in Fig. [3]jf). In 
the linear regime, the beam broadens and undergoes Zit- 
terbewegung |10[lll|. whereas in the nonlinear regime soli- 
ton formation is expected to take place with suppression 
of both beam broadening and Zitterbewegung. This is 
clearly shown in Fig. [3j which indicates the formation of 
the soliton during propagation with parameters as in Fig. 
[2] The evolution of field profiles \a2n\ and |aan-i| at even 
and odd waveguide positions is depicted in Fig. [3]ja) and 
|3][b), respectively. The evolution of the Fourier trans- 
form of the field a n of Fig. [3ja,b) along z is shown in Fig. 
(3jc). One can see that the strong component ai n in Fig. 



4 



[3^a) does not change much during propagation, whereas 
the weak component <Z2n-l m Fig- [31(b) is dramatically 
altered during propagation. As seen from Fig. ^b), at 
the beginning of the propagation the beam undergoes the 
Zitterbewegung. After reaching z ~ 70, the profile |ti2n— 1| 
becomes stable. Figure [3jd) shows the strong component 
|fl2n| °f the soliton profile with solid curves and the weak 
component |oan— 1| with dashed-dotted curves. As in Fig. 
[2]Jc,d) field profiles are taken at four values of propaga- 
tion distance z = (red curves); 50 (blue curves); 140 
(green curves); and 200 (black curves). One can also see 
that the strong component |a2 n | is stable, whereas the 
weak component first gets distorted (see blue and green 
curves) , but eventually the output curve (black color) re- 
laxes to the input curve (red color). Figure [3je) depicts 
the phase pattern of the field amplitudes across the ar- 
ray calculated at different z with corresponding colors 
as in Fig. J3^cl) . At the input (red curve) we have the 
phase difference equal to ir/2 between adjacent waveg- 
uides, but this phase pattern quickly transforms into the 
phase pattern of the soliton solution given by Eqs. 
i.e., S n = ...(p, p),{p + n, p + n),(p, />)... [see blue, green 
and black curves in Fig. [3][e)]. Therefore, here one can 
make a local conclusion: a beam with the intensity pro- 
files of soliton solution given by Eqs. ([6]), but with phase 
difference equal to tt/2 between adjacent waveguides will 
first undergo Zitterbewegung, but eventually its intensity 
profile and phase pattern will relax to the ones of the 
soliton solution given by Eqs. 

Dirac solitons - As mentioned in the introduction, 
BWAs have been used to mimic phenomena in both non- 
relativistic and relativistic quantum mechanics. To the 
best of our knowledge so far all these phenomena which 
have been simulated by BWAs are linear. In this section 
we will report on the simulation of nonlinear relativis- 
tic Dirac solitons in BWAs. As shown in [lOj [H] linear 
CMEs [Eqs. ([!])] for a beam with phase difference equal 
to 7r/2 can be converted into the linear one-dimensional 
relativistic Dirac equation (DE). Note that Eqs. ([!]) can 
be converted into the DE only for beams with special 
phase patterns; for instance, at normal beam incidence 
Eqs. ([I]) can not be converted into the DE. It turns out 
that with the soliton solution given by Eqs. (|6]), one 
can also successfully convert Eqs. ([l} into the nonlin- 
ear relativistic Dirac equation (NDE). Thus, one can use 
BWAs to mimic the relativistic Dirac solitons, and soli- 
ton solutions in BWAs given by Eqs. §6§ can be used 
to construct directly the Dirac soliton. Although the so- 
lution of Eqs. (|6| does not possess the phase difference 
equal to 7r/2 between adjacent waveguides [see Fig. |2jd)], 
the fact that it exhibits two wavenumbers k = ±ir/2 [see 
Fig. [2jb)] gives us some hope that the NDE can also be 
obtained in this case. Indeed, this is the case as shown 
below. In general, suppose that [a 2n (z), a 2n _i(z)] T = 
i 2n [g(2n, z),q(2n — l,z)] T , where the two functions g 
and q are smooth and their derivatives d n g and d n q ex- 
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FIG. 3: (Color online) (a,b) Propagation in the (n, z)-plane 
of the even and odd components of the beam with initial 
phase difference equal to it/2 between adjacent waveguides. 

(c) Fourier transform of field amplitudes in the (k, z)-plane. 

(d) Absolute values of the field amplitudes for intense com- 
ponent \a,2n\ with solid curves and weak component |<i2n-i| 
with dashed-dotted curves at four different values of z — 
(red curves); 50 (blue curves); 140 (green curves); and 200 
(black curves), (e) Phase pattern 5/n of field amplitudes for 
the same values of z as in (d). Colors of curves in (e) have 
the same meaning as in (d) . (f ) Scheme of the BWA structure 
for generating discrete solitons. Parameters: k =1; 7 =1; 
— 1.2; and no = 5. 

ist in the quasicontinuous limit [Eqs. ^ satisfy these 
requirements]. After setting ^(n) = ( — l) n a 2 „ and 
^2{n) = i{— l)"a 2 „_i, and following the standard ap- 
proach developed in [101 111 J we can introduce the continu- 
ous transverse coordinate £ •<-> n and the two-component 
spmor z) = (*i, * 2 ) T which satisfies the ID NDE: 

id z $ = -inad^ + a(3$ -jG, (7) 

where the nonlinear terms G = (|\E r i| 2 \l/i, |* 2 | 2 5' 2 ) T ; 
P = diag(l, —1) is the Pauli matrix a z ; and a is the 
Pauli matrix a x with diagonal elements equal to zero, 
but off-diagonal elements equal to unity. Note that Eq. 
is identical to the DE obtained in [TU1 EE] with the 
only difference that now we have the nonlinear term G 
in Eq. Q. Similar soliton solutions have been found for 
the NDE in Ref. [27], but with different and more com- 
plicated kind of nonlinearity, in the context of quantum 
field theory. Note that the nonlinearity that we have in 
Eq. violates Lorentz invariance [5S], and is similar to 
that of the Dirac equations in Bose-Einstein condensates 
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[T9] . Using the soliton solution given by Eqs. (JgJ and the 
above relation between a n and \& one can easily obtain 
the Dirac soliton solution of Eq. ([7| as follows: 



= sech(^-)e 

V V Tin ' 



iz(- 



-a) 



% ■> - 



sech(2^i)tanh(2£^i)6 



-cr) 



(8) 

The above solution is obtained for a > and 7 > 0. One 
can use the symmetry properties of Eqs. ([I]) to construct 
other Dirac soliton solutions of Eq. ([7]), with different 
sign combinations between a and 7. The expressions 
given by Eq. ^ give the main result of this Letter, and 
the only physically realizable way that we are aware of 
to produce and observe Dirac solitons with a table-top 
experiment. In future investigations we are planning to 
carefully study the dynamics and the stability of Dirac 
solitons in BWAs, on which we will report in a separate 
publication. 

Conclusions - In this Letter we have provided an- 
alytical expressions for the non-moving gap solitons in 
BWAs and shown their connection to Dirac solitons in a 
nonlinear extension of the relativistic ID Dirac equation 
describing the dynamics of a freely moving relativistic 
particle. Our results suggest that BWAs can be used 
as a classical simulator to investigate relativistic Dirac 
solitons, enabling to realize an experimentally accessi- 
ble model system of quantum nonlinearities that have 
been so far a subject of speculation in the foundation 
of quantum field theories. The analysis of analogue of 
quantum field theory effects as those ones described in 
this Letter is applicable to virtually any nonlinear dis- 
crete periodic system supporting solitons, either classical 
or quantum, therefore making our results very general 
and of relevance to different systems beyond optics, such 
as ultracold atoms in optical lattices and trapped ions 
where analogs of linear relativistic effects, such as Zitter- 
bewegung, have been studied and observed [28 30 
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